Treed Gaussian processes for animal movement modeling

Abstract Wildlife telemetry data may be used to answer a diverse range of questions relevant to wildlife ecology and management. One challenge to modeling telemetry data is that animal movement often varies greatly in pattern over time, and current continuous‐time modeling approaches to handle such nonstationarity require bespoke and often complex models that may pose barriers to practitioner implementation. We demonstrate a novel application of treed Gaussian process (TGP) modeling, a Bayesian machine learning approach that automatically captures the nonstationarity and abrupt transitions present in animal movement. The machine learning formulation of TGPs enables modeling to be nearly automated, while their Bayesian formulation allows for the derivation of movement descriptors with associated uncertainty measures. We demonstrate the use of an existing R package to implement TGPs using the familiar Markov chain Monte Carlo algorithm. We then use estimated movement trajectories to derive movement descriptors that can be compared across individuals and populations. We applied the TGP model to a case study of lesser prairie‐chickens (Tympanuchus pallidicinctus) to demonstrate the benefits of TGP modeling and compared distance traveled and residence times across lesser prairie‐chicken individuals and populations. For broad usability, we outline all steps necessary for practitioners to specify relevant movement descriptors (e.g., turn angles, speed, contact points) and apply TGP modeling and trajectory comparison to their own telemetry datasets. Combining the predictive power of machine learning and the statistical inference of Bayesian methods to model movement trajectories allows for the estimation of statistically comparable movement descriptors from telemetry studies. Our use of an accessible R package allows practitioners to model trajectories and estimate movement descriptors, facilitating the use of telemetry data to answer applied management questions.

To work through this appendix, we use 100 consecutive GPS points from one lesser prairie-chicken.After using this example to follow the code it can be replaced with your own data and research questions/derived quantities.
We rely on Robert Gramacy's tgp package (version 2.4) to fit movement models.

#install.package("tgp") library(tgp)
Step 1) Read in data The dataframe we use requires X (GPS easting or longitude), Y (GPS northing or latitude), and t.obs (GPS recorded time) columns, which may require some manipulation of your GPS data.Data should already be cleaned and checked for location error outliers, as described in Appendix S3 section 3.
Modeling is performed either using lat-long or UTM coordinates; predicted trajectories will be outputted in whichever form the coordinates are read in as.This means the user may decide which coordinate form to use.If modeling in lat-long, trajectories will have to be converted to UTM to perform derived quantity computations such as distance between points.However, UTM may pose problems for large datasets that span multiple zones.The lesser prairie-chicken dataset utilizes UTM.

Add t.obs column:
t.obs provides a standardized time measure for the model to intake.Units will depend on your time scale, but counting in hours is common.Hour "0" can relate to any time within or outside the dataset, and won't affect the modeling, as you'll transform back to the natural time scale for your results.We used t.obs = hours since first recorded GPS point.#convert to hours since start individual$t.obs<-(as.numeric(individual$timestampmin(individual$timestamp)))/3600 Step 2) Fit TGP Model

Fit the model to the X (easting) and Y (northing) coordinates separately
Bayesian treed Gaussian process with jumps to the limiting linear model (within the tgp package) Note that for large datasets and without any parallelization this may be a time consuming step.
Inclusion of jumps to the limiting linear model in the TGP model allows the model to fit a linear model in partitions where the data do not require a Gaussian process, see Gramacy & Lee (2008)  Step 3) Predict paths

Set MCMC hyperparameters
Prediction under the treed Gaussian Process requires MCMC sampling.In the tgp package this requires the tuning parameter of BTE = c(burn in, total samples, thinning value).Setting of these MCMC sampling parameters is explained in Appendix S3 section 2, and we recommend an understanding of MCMC sampling at the level of Hobbs and Hooten (2015) or Hooten and Hefley (2019).
We recommend prototyping your data with small BTE settings first, for example c(400,500,2

Use fitted model to predict
This chunk does not require editing.It will sample from the posterior predictive distribution of locations at each ∆t time point.
Note that for large datasets and without any parallelization this may be a time consuming step and should be prototyped on a small dataset and small BTE first.After all_paths is outputted, this would be a good time to export all_paths and save it on your device if your dataset is large.
Though it may be tempting to elimininate this storage required for the all_paths data frame by predicting the paths and computing the derived quantities simultaneously, saving the full MCMC sampled distribution of fine-resolution trajectories allows the flexibility to change scale or derived quantities as research questions change or are added.Because predicting these full trajectories is the most time consuming step, we recommend saving the full predicted trajectory sample outside of R for future use.
To check the effective sample size of your MCMC samples, run the following code after you've sampled from the posterior predictive distribution.A large effective sample size (close to the MCMC sample size) implies good mixing of the sample and low autocorrelation.Low effective sample size suggests autocorrelation that may require increased thinning (increased E in BTE settings).

## var1 ## 1000
To examine a trace plot of your MCMC samples, run the following after you've sampled from the posterior predictive distribution.

Plot predicted trajectory
We can visualize the predicted trajectory in both one and two dimensions.Black squares plot recorded GPS data points and blue circles plot the mean (expected value) of the MCMC samples from the posterior distribution of predicted locations at the set ∆t inteval.Plotting code can be found in the Appendix S2.Rmd file One dimension, including 95% credible intervals (gray lines) from the full MCMC sample: Step 4) Compute derived quantity at ∆t resolution First we use the all_paths table to compute the derived quantity at the scale of ∆t (e.g., hourly).This is the scale that all_paths is already at.Beginning at this resolution allows us to then transform to any larger desired scale (e.g., daily, weekly).
It may be helpful to follow along with the displacement example in Appendix S4 as well, and we've noted how each code step aligns with each step in Appendix S4.

Your derived quantity
You will have to edit this section depending on your derived quantity.As the options for derived quantities are infinite, this may take some creativity.Refer to Table 1 of the manuscript for basic derived quantity functions that can be built off of.
You will have to write your derived quantity function here: The displacement function used in section 4.2 of the manuscript is provided as an example.We now have a data frame of samples of the derived quantity (in this case hourly displacement) at each ∆t point, with d1-d1000 indexing the 1000 samples.

Step 5) Output inference at desired scale
The above 252,000 samples of different hourly displacements may not be at the temporal scale of interest, so we can transform it into our desired value: average daily total distance traveled over the study period.
Though it may be tempting to just average all values in the table, we must maintain alignment of MCMC samples (see Chapter 8.3 in Hobbs & Hooten, 2015).See Appendix S4 for further walking through of this example.

These steps will vary based on your transformations and desired ending derived quantity
In our example, if we want to transform to daily total distance traveled averaged across the season, we first take daily sums within each of the MCMC samples (columns) separately.
[arrow 2 in Appendix S4] We used the dplyr and lubridate packages to make this easier.This gives a distribution (1000 samples) of daily total distance traveled.
To get 1000 samples of daily distance traveled averaged across the season, we take column (sample) averages.For inference at a different resolution, the above chunks must be edited.

Step 6) Repeat for individuals and study periods
We leave it to the practitioner to repeat this across individuals and chosen study periods, and to store the all_paths tables, and their desired derived quantity point estimates and variance estimates.

Speed up options
There are multiple options to increase run speeds.
Option A uses simple parallelization to fit both X and Y models simultaneously on your machine by using separate cores simultaneously.This can expect a modest ~50% speed up for only the model fitting step.
Option B offers the largest speed up (up to 10x), but requires slight reconfiguration of the computer and will depend slightly on the computer, making it more complicated.
Option C offers 2x or 3x speed up and is the most complicated to install, and requires reconfiguring the tgp package.
If speed is an important goal, combining option B and option C will have the best outcome.

Option A)
Simple parallelization of fitting the X and Y models.(This simple parallelization cannot be applied to prediction.)

library(parallel)
Replace fitModels chunk with: Option B is to use a more efficient matrix algebra library within R. This will greatly speed up computations done within both stages of tgp.These changes must be made on the computer outside of R, but once this change is made then tgp fitting and prediction will automatically run faster in R every time, without any changes to code.This can be done in the Accelerate framework on a mac, following installation instructions provided by Berkley statistics.
For windows, Intel's Math Kernel Library (MKL) can be used.
More details on Accelerate and MKL can be found in Robert Gramacy's Surrogates (2020).

Option C)
Option C requires uninstalling and reinstalling tgp with Pthreads.This will make tgp prediction automatically utilize parallelization.Details for this and reinstillation can be found in the tgp documentation appendix C.2.
Two dimensions, including 20 sampled trajectories in gray to visually represent the MCMC sample: